The present invention relates to a method of controlling the quality of seismic data migration by using the generalized Radon transform (GRT).
There exist numerous techniques that are used by geophysicists to find out about the geological layers constituting the subsoil. One of these techniques, known as a xe2x80x9cwell seismic surveyingxe2x80x9d is shown in FIG. 1 and consists in emitting sound signals at the surface of the ground (source S) and in using sensors Ci located in a well 10 to record the seismic waves as they are transmitted after reflection and/or refraction on the boundaries between the various geological layers. After the waves have been recorded by the sensors, they are filtered to determine which waves were reflected by reflectors situated beneath the sensors. FIG. 2A shows the signals as recorded prior to filtering and FIG. 2B shows the corresponding signals after filtering. A migration method applied to the filtered signals then makes it possible to obtain a migrated image which can be thought of as a vertical section plane through the subsoil containing the source(s) and the sensors in the well. As shown in FIG. 3. the migrated image is constituted by a succession of vertical recordings referred to as xe2x80x9cmigrated seismic tracesxe2x80x9d constituted by signals of greater or lesser amplitude made up of positive and negative arches. The signals relating to the same reflector correlate from one trace to the next, so that the eye can see delineations which correspond to boundaries between geological layers.
The document xe2x80x9cThe generalized Radon transform, a breakthrough in seismic migrationxe2x80x9d, by M. Oristaglio, G Beylkin. and D. Miller (Seismics, The Technical Review. Volume 35, Number 3, 1987. pages 20 to 27) describes an example of the migration method which makes use of the generalized Radon transform (GRT). The GRT migration method constructs the migrated image point-by-point. As shown in FIG. 4, each image point I1 or I2 is considered as a point of velocity change in the medium that is capable of reflecting energy coming from the source S. The energy diffused by said point propagates in all directions and is finally detected by the sensors. The travel time of the seismic wave is calculated by dividing the sum of the distance ls from the source to said point plus the distance lc from the point to the sensor by the assumed propagation speed of the wave in the subsoil. This propagation speed is referred to as the velocity of the medium.
For each particular reflector point in the space of the migrated image it is possible to obtain a travel time curve in the data space: CI1, CI2; by calculating the seismic wave arrival times at each sensor as shown in FIG. 5.
Inversely, each point P in the data space corresponds to the amplitude of the signal measured by a sensor at a given instant, which measured signal is due to a plurality of waves reflected by different reflector points all reaching the sensor at the same time. Thus, every point of the migrated image whose travel time curve passes through a given point P contributes to the amplitude of the signal at that point. The set of these reflecting points constitutes an isochronous curve which. as shown in FIG. 4, is elliptical in shape, with the source S and the sensor C being at the foci, assuming that the velocity in the subsoil is constant.
In application of the generalized Radon transform, each point of the migrated image is obtained by summing the contributions of data along the corresponding travel time curve, said contributions being weighted to take account of physical phenomena relating to wave propagation, for example loss of amplitude proportional to path length. To calculate path curves by ray tracing, an initial speed model is used which integrates a priori knowledge about the subsoil, such as the model shown in FIG. 6.
The GRT migration algorithm thus has the following steps:
For each data trace T, for each point (x, z) of the migrated image:
i) the travel time ts from the soure S to the point (x, z) is calculated;
ii) the travel time tc from the sensor C to the point (x, z) is calculated;
iii) this gives the source-to-sensor travel time: tsc=ts+tc;
iv) the sample u of the trace T at time tsc is found;
v) the-weighting factor w relating to the source S, to the sensor C, and the point (x, z) is found;
vi) the value of an image table (x, z) at said point is calculated;
Image(x, z)=Image (x, z)+uxc3x97w.
A variant of this GRT algorithm consists in a first step in calculating the table of travel time parameters t, the path length 1, and the direction of the ray at its ends d and a, and in a second step in performing the migration calculation. Thus, as shown in FIG. 7, for each source and for each sensor C, the travel times t(P, Ixy) to go from said source or said sensor to all of the migrated image points Ixy are calculated. The path lengths l(P, Ixy) and the orientations of the ray at its ends a(P, Ixy) and d (P, Ixy) are calculated. In the second step. migration uses the table calculated during the first step.
For example, with filtered seismic signals as shown in FIG. 8, GRT migration enables a migrated image to be obtained of the kind shown in FIG. 9.
Before using the results as shown in FIG. 9, the operator must verify that they are well-founded. It is necessary to perform quality control thereon. Traditionally. this has been based on more or less subjective criteria such as how the results correspond with a priori knowledge about the subsoil, or the continuity of the delineations representing the reflectors. It is thus very difficult to perform such quality control. particularly since the representation of the data before migration (FIG. 8) is very different from the representation of the data in the migrated image (FIG. 9).
A method described in xe2x80x9cAutomatic association of kinematic information to pre-stack imagesxe2x80x9d, by S. Geoltrain and E. Chovet (61st Annual International Meeting, Sot. Expl. Geophys., Expanded Abstracts, 1991, pages 890 to 892), for Kirchhoff migration. makes it possible to associate each point of the migrated image with a point in data space. Nevertheless, that method does not enable each point in the migrated image to be associated with all of the contributing points along the travel time curve. Nor does it make it possible to find the isochronous curve associated with a point in data space.
An object of the present invention is to solve those problems by proposing a quality control method for GRT migration of seismic data by interactive identification of events between the image obtained by migration and the seismic data prior to migration.
The present invention provides a method of controlling the quality of seismic data that has been migrated using the generalized Radon transform, the method serving to pass between a seismic data space and a migrated image space, in which method a parameter table is calculated giving, for a seismic wave going between a point of the image and a source or a sensor, its path length, its travel time, and the angles of incidence of the wave at the beginning and at the end of the path, the method being characterized in that correspondence is established between at least one zone of a first one of said two spaces and at least one zone of the second space by using said parameter table to fill in a correspondence table QCimage.
The parameter table is usually obtained by ray tracing, but it can also obtained by any other means making it possible to determine, for a seismic wave going between an image point and a source or a sensor, its path length, its travel time, and the angles of incidence of the wave at the beginning and end of the path.
The method of the invention makes it possible firstly to point to a zone in seismic data space and to project it into migrated image space by local migration relating to the pointed-to zone (xe2x80x9cforward projectionxe2x80x9d). The method also makes it possible to point to a zone in migrated image space and to project it into seismic data space by calculating the seismic data field that contributes to the zone pointed to in the migrated image (xe2x80x9cbackward projectionxe2x80x9d).
In the first embodiment (forward projection) the point (x, z) of the migrated image is put into correspondence with the corresponding zone in data space if the value QCimage(x, z) of the correspondence table is greater than a given threshold, e.g. zero.
Advantageously, the zone pointed to in data space is a set of intervals Zi each characterized by a trace T. and the start and finish td and tf of the corresponding interval. For each interval Zi and for each point (x, z) of the migrated image, the following steps are performed:
i) for the ray to the sensor C relating to said interval Zi at the point (x, z). the following parameters are read from the parameter table: travel time tc, path length lc, and orientation of the ray at its ends dc and ac;
ii) the source-to-sensor travel time tsc tc+ts is calculated;
iii) if tdxe2x89xa6tscxe2x89xa6tf, then the following weighting factor w(ls, ds, as, lc, dc, ac) is calculated;
iv) the value QCimage(x, z)=QCimage(x, y)+w is calculated;
v) loops are performed for the points and the intervals so as to fill the correspondence table QCimage and
vi) said table QCimage is used for establishing correspondence between zones.
In the second embodiment (backward projection), the point corresponding to a trace T at source-sensor travel time tsc of the data is put into correspondence with the corresponding zone in the migrated image if the value QCimage (T, tsc) in the correspondence table is greater than a given threshold, e.g. zero. Advantageously, the zone pointed to in the migrated image is a set of intervals Zi characterized by the position of the trace xz, and the start and finish zd and zf of the corresponding interval. For each trace T of input data, for each interval Zi, and for each point (x, z) of the interval Zi, the following steps are performed:
i) for the ray from the source S relating to the trace T at the point (x, z), the following parameters are read from the parameter table: travel time ts, path length ls, and orientation of the ray at its ends ds and as;
ii) for the ray to the sensor C relating to the trace T at the point (x, z), the following parameters are read from the parameter table: travel time tc, path length lc, and orientation of the ray at its ends dc and ac;
iii) the source-to-sensor travel time tsc=tc+ts is calculated;
iv) the weighting factor w (ls, ds, as, ls, dc, ac) is calculated;
v) the value of QCimage(T, tsc)=QCimage(T, tsc)+w is calculated;
vi) looping is performed for the points, the intervals,and the traces to fill the correspondence table QCimage; and
vii) said table is used to establishes correspondence between zones.
Advantageously, in said method a parameter whose value is proportional to the value of QCimage is assigned to the corresponding zones in the two spaces. By way of example, this parameter can be a semitransparent color, and its value can be the density of the color.
The method of the invention is usable in numerous fields,in particular in well seismic surveying, and also in xe2x80x9cpre-stackxe2x80x9d surface seismic surveying, with the traces coming from various sources and sensors not being added to one another in this case, thus making it possible to one another in this case, thus making it possible to apply the GRT migration method.